Elastic critical behaviour in a 3d model for polymer gels 
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, O, . The elastic response in polymeric gels is studied by means of a percolation 

cn . dynamic model. By numerical simulations the fluctuations in the gyration 

^^ , radius and in the center of mass motion of the percolating cluster are de- 

j~^ , termined. Their scaling behaviour at the gelation threshold gives a critical 

o. 

exponent for the elastic modulus / ^ 2.5 it 0.1 in agreement with the predic- 



tion f = dv. 
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I. INTRODUCTION 

The elastic response in polymer gels is due to the presence of the polymeric network 
formed in the gelation process: this random structure makes the system able to respond to 
external stresses. Due to the large possibility of conformational changes which characterizes 
this macromolecular phase, there is an entropic term contributing to the elastic properties 
of the system. This picture is typical of gelling systems and may fit a wide range of very 
different materials. Then because of the different role played by the entropic term and the 
particular features of the polymeric networks, a rich phenomenology is produced. 
The viscoelastic behaviour of gelling materials and the mechanisms at its origin certainly 
represent a central issue in soft matter physics and have a high relevance to a wide range 
of applications from food processing to materials science. On the other hand these mate- 
rials show many typical features of complex systems and their study is connected to some 
fundamental problems in soft matter physics, from the entropic elastic behaviour, which we 
are interested in here, to glassy dynamics. As a consequence these systems are intensively 
investigated both experimentally and by means of statistical mechanics models, bringing up 
a lively debate. 

In particular here we study the critical behaviour of the elastic response of the system as 
the gelation process takes place producing the polymeric network: in a gelling system the 
elastic modulus starts growing at the gelation threshold with a power law behaviour, usually 
expressed in terms of the polymerization degree. In the experiments performed on different 
gelling materials the value of the critical exponent / describing the critical behaviour of 
the elastic modulus appears to be close to either / ~ 2, or / ~ 3, or else / ~ 4. The 
value / ~ 2 is actually observed in experiments on agarose gel [1], gelatin gels [2], some 
silica gels (TEOS) [3]. Within statistical mechanics models it corresponds to the prediction 
based on the de Gennes' analogy between the elasticity of a percolating random network 
of Hookean springs and the conductivity in a random percolating network of resistors [4]. 
The critical exponent / ~ 3. has been observed in diisocianate/triol gel [7], in epoxy resins 



[8], TMOS silica gels [9], polyesters gels [10]. It is very close to the value / = du, where 
d is the dimensionality of the system and u is the critical exponent of the connectedness 
length C,, predicted in ref. [11] which in 3d gives / ~ 2.64 (in mean field both du and the 
random resistor network exponent are equal to 3). Following the argument in ref. [11] this 
exponent would then correspond to the case of a dominating entropic term, evaluated by 
means of scaling arguments for the percolating network: the elastic energy of the system is 
due to the entropic elasticity of the macrolink whose length is of the order of the connect- 
edness correlation length. Finally the value / ~ 4 can be linked to the prediction of the 
bond-bending model f = du + 1 [12,13] or alternatively f = t + 2v [14,15], which in 3d give 
/ ~ 3.7. This would imply that within the elastic energy describing the system there is a 
bending term playing a relevant role in the elastic response of the network. This critical 
behaviour is typical of some colloidal gels [16]. The clustering of the experimental values 
for the elasticity critical exponent in gels around discrete values suggests the possibility of 
correspondently individuating different universality classes, which should be characterized 
by some intrinsic features of the networks formed in the materials. 

Recent numerical studies via molecular dynamics [5] and Monte Carlo simulations [6] of 
percolating networks of tethered particles with no hard core interactions have shown that 
the shear modulus critical exponent is ~ 1.3 in rf = 2 and ~ 2. in (i = 3. Monte Carlo 
simulations of two and three-dimensional percolating networks of tethered particles with 
hard core repulsion [6] find consistent results for the shear modulus critical exponent. These 
results agree with the de Gennes prediction and with some experimental results. 
Within the numerical studies we have approached the study of this problem introducing a 
percolation dynamic model, and directly investigating the dynamic viscoelastic properties 
as the percolation transition takes place. Our model introduces in the percolation model 
the bond-fluctuation dynamics, which takes into account the conformational changes of the 
polymer molecules and the excluded volume interactions. This model has been translated in 
a lattice algorithm and studied via numerical simulations on hypercubic lattices. Actually it 
presents many fundamental features of the gelation phenomenology and has already allowed 



to study the critical behaviour of the viscoelastic properties and the relaxation process in 
gelling systems in the sol phase [19]. Numerical simulations of the model in d = 2 [20] have 
shown that the elasticity critical exponent is ~ 2.7, a value which agrees with the prediction 
f = du in [1 1] . The study of the model in rf = 3 is of great interest to check the elasticity 
critical behaviour and to eventually compare the findings with the experimental results. We 
here present the results of large scale 3d simulations: the data show that / ~ 2.5 ± 0.1, 
again in agreement with the prediction / = du, coherently with the findings in d = 2. 
In the next section we present the model and the details of the numerical simulations, in sec- 
tion III the elastic response of the percolating cluster is discussed and a scaling behaviour is 
obtained; in section IV the results of the numerical simulations are presented and discussed; 
section V contains concluding remarks. 

II. MODEL AND NUMERICAL SIMULATIONS 

We consider a solution of tetrafunctional monomers of concentration p. The monomers 
interact via excluded volume interactions, i.e. a monomer occupies a lattice elementary 
cell and two occupied cells cannot have common sites. Nearest neighbours or next-nearest 
neighbours are instantaneously linked by a permanent bond with probability pi,. In terms 
of these two parameters the percolation line can be determined via the critical behaviour of 
the percolation properties, individuating the sol and the gel phase in the system. Actually 
in the simulation we fix pf, = 1 and study the system varying p [19-21]. The percolation 
quantities critical exponents agree with the random percolation predictions [18] (e.g. 7 ~ 
1.8 ± 0.05 and u ~ 0.89 ± 0.01 in 3d [19]). The monomers free or linked in clusters diffuse 
via random and local movements on the lattice according to the bond-fluctuation dynamics 
[22], which is ruled by the possibility of varying the bond lengths within a set of values 
determined by the excluded volume interactions and the SAW condition. This produces a 
high number of different bond vectors and we consider the case of permanent bonds, which 
corresponds to the strong gelation process. In 3d the allowed bond lengths on the cubic 



lattice are / = 2, ^/5, ^/6, 3, v^TO and in Fig.l an example of different allowed configurations 

for a polymer molecule is shown. 

We here present the results of extended numerical simulations of the model in the gel phase 

to study the elastic response in the system. It is worth noticing that there is no elastic 

potential energy for the bond vectors and then the elastic behaviour is purely entropic 

implying that our study is performed at finite temperature. 

The data presented here refer to lattice sizes L ranging from 12 to 32, and have been averaged 

over 30 different realizations. The simulations have been performed on the CRAY — T3E 

system at CINECA taking more than 30000 hours/node of CPU time. 

III. ELASTIC RESPONSE IN THE GEL PHASE 

We study the elastic response in the gel phase in terms of the macroscopic elastic constant 
of the system K, which is experimentally defined as the ratio between an applied external 
force and the deformation. In a simple elongation experiment if Iq is the undeformed length 
and 6 = {I — /q) is the deformation in the system, within the linear response approximation 
the elastic free energy F ~ K6^. In terms of the Young elastic modulus E the free energy 
per unit volume is F/V ~ E6'^/ll. Then K ~ EV/ll and for a cube of size L, K r^ EL'^~'^, 
expressing the fact that the elastic modulus has the dimensions of an energy per unit volume 
and is an intensive quantity whereas K depends on the system size L. 
In the gel, since E vanishes at p^ as ~ .^^-^ (where / = f/iy) one has 

K ~ L'^-^C^ (1) 

Following eq.(l) the macroscopic elastic constant presents the corresponding scaling be- 
haviours as function of the system size L and of the distance from the percolation threshold 

(P-Pc)- 



P > Pc K r^ L'^ "^ 

p = p, K^ L-' (2) 

fixed L K r^ (p — pc)^ 

where z = f — {d — 2). 

An alternative way to obtain these scahng relations is to consider the percolating cluster as 
a network of nodes connected by macrolinks of linear size C, [23]. Each macrolink can be 
considered as a spring with an effective elastic constant. At the percolation threshold there 
is only one macrolink spanning the system and its effective elastic constant coincides with 
the system macroscopic elastic constant K. 

In order to evaluate K we notice that for a spring of elastic constant i^ in a thermal bath 
at temperature T, the mean fluctuation in the energy U is (Af/) = ^K{x'^), where x is the 
spring elongation. From the energy equipartition K{x'^) = ksT, so that at the equilibrium 
the elastic constant K is related to the the fluctuations in the spring length [24-27]. This 
result can be more generally obtained by means of the Fokker-Plank equation for the prob- 
ability distribution of the spring elongation [24]. 

Therefore the macroscopic elastic constant of the gel phase is related to the fluctuations of 
the linear size of the infinite cluster, i.e. the squared fluctuations in the gyration radius of 
the percolating cluster (Ai?2) = ((/?^ _ {Rg)^), K ~ l/(Ai?2). 

An alternative way to calculate the macroscopic elastic constant in the gel phase by means 
of the fluctuations in the unperturbed system is to consider the center of mass of the perco- 
lating cluster as a brownian particle, subject to a restoring force responsible for the elastic 
behaviour of the system. The restoring force introduces a limitation on the diffusive pro- 
cess of a brownian particle. Using the same argument based on the energy equipartition, 
the asymptotic equilibrium value A of its displacement fluctuations (Ai?^(t)) is inversely 
proportional to the elastic constant, i.e. A ~ 1/ K. 



IV. RESULTS 

In order to numerically study the elastic response we have used both the approaches 
mentioned before namely we have calculated the average fluctuation of the gyration radius 
of the percolating cluster (Ai?^ and the asymptotic value A of the mean square displacement 
of its center of mass {AR^{t)) [28]. 

In the flrst approach the average fluctuation in the gyration radius (Ai?^ of the percolating 
cluster has been computed at the percolation threshold in system of different size L with 
periodic boundary conditions. In Fig.2 (Ai?p at Pc is shown in a log-log plot as function 
of the lattice size L. In the considered range the data are well fltted by a power law 
behaviour according to eq.(2), giving a critical exponent 5 ~ 1.9 ± 0.1, i.e. z ~ 1.7 ± 0.1. 
As z = f - {d - 2), f = z + {d - 2)u and this resuh gives / ~ 2.6 ± 0.1. 
In the second approach we have calculated (Ai?^(t)) at different steps of the gelation process, 
i.e. as p grows above the percolation threshold Pc ~ 0.718±0.005 [19], and at Pc for different 
lattice sizes. In these simulations hard- wall boundary conditions have been used [29]. In 
the gel at the gelation transition the percolating cluster is a quite loose network, the center 
of mass is rather free and the elastic response is weak. As the gelation process goes on, 
the network tightens becoming more rigid, the elastic constant of the system increases and 
the center of mass motion is progressively constrained. This would be then the physical 
mechanism producing the critical behaviour of the elastic response for a critical gel. In 
agreement with our picture (Ai?^(t)) grows with time up to a limiting plateau value A as 
it is shown in Fig. 3. This quantity is inversely proportional to the elastic constant of the 
system A ~ 1/K and increases as the percolation threshold is approached from above. 
In Figures 4,5 and 6 the scaling behaviour obtained for A is presented. 
In Fig. 4 A{L,p) for p >> pc {p = 0.85) is shown in a double logarithimic plot: being far 
from the percolation threshold the system is reasonably homogeneous and in fact the scaling 
behaviour as A ~ /^-o.99±o.i jg observed, in agreement with the behaviour K ~ L'^^'^. 
In Fig. 5 A[L,p = p^) is shown in a double logaritmic plot as a function of the lattice size L: 
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the data exhibit a behaviour A ~ L^ with z ~ 2.0±0.1. This result again gives / ~ 2.6±0.1. 

Finally in Fig.6 A(p — pc) is plotted for the lattice size L = 32: we fit the data with a power 

law behaviour (eq.(2)) and obtain a critical exponent / ~ 2.5 ± 0.1. 

It is straightforward to notice that all these numerical results can be coherently interpreted 

in terms of the scaling relations obtained for K. 

The value of the critical exponent / ~ 2.6 is in good agreement with the prediction f = du 

of ref. [11], therefore supporting the picture there proposed, and consistent with the value 

obtained in the 2d study of the model [20]. 

Due to the limited extension of the critical parameter {p — pc) and of L here investigated 

the eventual occurrence of a crossover to a different exponent cannot be excluded. 

V. CONCLUSIONS 

The numerical results of Figures 2,5 and 6 show that z = 2., coherently agreeing with 
the prediction / = du. On the whole they support the scaling picture we propose and the 
argument of ref. [11]. They also agree with some experimental results [7]- [10]. This result 
has been obtained via two independent calculations giving consistent numerical values and 
is also consistent with the value previously obtained in the 2d study. 

On the other hand the recent numerical works on entropic elastic models of refs. [5,6] find 
a good agreement with the de Gennes prediction. These results, together with the experi- 
mental data, seem to suggest the possibility that there are two distinct universality classes 
characterized one by an exponent f = du and another by the electrical analogy exponent 
f = t, which in 3d are respectively ~ 2.64 and ~ 2.. However since the models in the 
different numerical studies are rather similar the possibility of a crossover between different 
dynamic regimes, as it is observed in some experiments [3], cannot be completely excluded. 
In both cases these results give a hint for the interpretation of the experimental data and 
indicate the aspects which should be further investigated. 
This work has been partially supported by the European TMR Network- Fractals under con- 
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FIGURES 
FIG. 1. An example of time evolution of a cluster formed by four monomers according to the 

bond-fluctuation dynamics: in a, starting from the upper central bond and clockwise, the bond 

lengths are I = v5, 3, 3, 2; in b the upper left monomer has moved forward and / = 2,3, 3, v5; 

having moved right the other left monomer in b one has c with I = 2,3, \/^, -v/6; moving right the 

front monomer in c the d configuration is obtained with 1 = 2, vTO, \/6, \/6. 

FIG. 2. Log-log plot of the fluctuation of the percolating cluster gyration radius (Ai?p as a 
function of the lattice size L at pc'- from the fit of the data the critical exponent z ~ 1.9 it 0.1 is 
obtained. Here and in the following figures the lengths are expressed in units of lattice spacing. 

FIG. 3. The mean square displacement of the center of mass of the percolating cluster 
(Ai?^(t)) as function of time for different value of the monomer concentration p: from below 
p = 0.8,0.77,0.76,0,75. The dashed lines correspond to the asymptotic plateau values A which 
grows approaching pc- The data refer to a lattice size L = 32. The unit time here is the MonteCarlo 
step per particle. 

FIG. 4. Log-log plot of A{L,p) for p >> pc (the data refer to p = 0.85) as function of L: the 
data show a scaling behaviour ~ ]^y'j;^0.99±0.i_ 

FIG. 5. Log- log plot of the plateau values A(L,p = Pc = 0.718) as function of L: the data are 
fitted by a power law giving the critical exponent z ~ 2.0 it 0.1. 

FIG. 6. The plateau values A(L = 32, {p — pc)) in a double logarithmic plot as function of 
(p — Pc)'- the best fit of the data close to the percolation threshold gives the critical exponent 
/~ 2.5 ±0.1. 
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